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ABSTRACT 

A finite element method is used to solve the equations of motion for 2- 
and 3-D fluid flow. The time-dependent equations are solved explicitly 
using quadrilateral (2-D) and hexahedral (3-D) elements, mass lumping, 
and reduced integration (when appropriate) . A Petrov-Galerkin technique 
is applied to the advection terms. The method requires a minimum of 
computational storage, executes quickly, and is scalable for execution 
on computer systems ranging from PCs to supercomputers. 
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W weighting function 

a Petrov-Galerkin parameter 

(3 Petrov weighting factor 

y ratio of specific heats 

At time step 

k thermal diffusivity 

p density 

a xx x-coordinate viscous stress 

a x x-y viscous shear 

o xz x-z viscous shear 

0 yy y-coordinate viscous stress 

0 y 2 y-z viscous shear 

a zz z-coordinate viscous stress 

D area (2-D) ; volume (3-D) 

Subscripts 

i,j,k local node index 
x x direction 

y y direction 

z z direction 

03 free stream 

column vector 

Superscripts 

T transpose 

-1 inverse 

n previous time value 

n+1 new time value 

x x direction 

y y direction 

z z direction 

u trial value 

1 1 average value 

INTRODUCTION 

Numerical modeling of multidimensional fluid flow is a difficult, 
computationally intensive effort. When the physical domain is simple, 
finite difference methods (EDM) are typically used as a matter of conve- 
nience and historical preference. If the domain is irregular, boundary 
fitted coordinates are usually employed [1]; however, such techniques 
are not easy to use, particularly when generating 3-D grids, and require 
the transformation of the governing equations. Finite element methods 
(FEM) , while mathematically more robust and convenient to use in complex 
geometries, can tax the limits of the largest computers. In an effort to 
alleviate these constraints, a modified FEM, which is competitive with 
the FDM in storage and speed, is used to solve a set of example problems 
in 2- and 3-D for both incompressible and compressible flows. The FEM is 
based on the method discussed in Pepper and Singer [2] and Pepper and 
Humphrey [3] . 
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The modified FEM employs simple, well known alterations to the 
conventional FEM: mass lumping, reduced integration, and application of 
Petrov-Galerkin weighting. Lumping the time derivative term into one 
diagonal allows simple explicit time integration to be used and the mass 
matrix to be easily inverted. One point quadrature, as opposed to 2x2 
(bilinear quadrilateral) or 2x2x2 (hexahedral) Gauss points, yields 
integral values obtained at the element centroid. This procedure permits 
explicit evaluation of the integral terms on a local level, i.e., per 
element without regard for bandwidth. The application of reduced quadra- 
ture has been found to be very effective when the element is not too 
distorted [4]; however, care must be exercised in generating the initial 
mesh to insure aspect ratios and distortions are minimal. 

Velocities within the advection terms are averaged over an element 
and a Petrov-Galerkin formulation applied [5] . Application of the 
Petrov-Galerkin technique helps maintain stability but does not overly 
dampen the solution, compared to the more commonly utilized upwinding 
techniques and anti-dispersive schemes applied in FDMs [6], 

Performance of the algorithm on current scalar/vector machines is 
good [2]; many multidimensional problems can be run on low-end comput- 
ers. The performance of the FEM is significantly improved on a parallel 
computer. Since the method is explicit in time, vectorization and 
parallelization are relatively easy. A version of the method has been 
written and optimized for execution on an Alliant VFX/40 mini- 
supercomputer with 4 processors. 

GOVERNING EQUATIONS 

The mathematical relations which describe the transport of fluid 
motion and energy are the equations for conservation of mass, momentum, 
and energy. The governing equations for three-dimensional flow can be 
written in non-dimensional form as [6] 

Conservation of Mass 

<3p dpu dpv dpw „ 
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dt dx dy dz 

Conservation of Momentum 
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(5) 


Conservation of Energy 
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where p is density, u is horizontal velocity, v is lateral velocity, w 
is vertical velocity, p is pressure, is total energy, q K , q y , and q z 
are the gradient flux terms for heat conduction, and o xx ,o xy ,o xz , etc., 
are the normal and tangential viscous stress terms. 
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( 6 ) 


dv 
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where R e = p„U a L/[i m . For compressible flow, pressure and temperature are 
obtained from the equation of state 

p = (y- l)pe 


7 = 


Y Mlp 


(9) 

( 10 ) 

is the free stream Mach num- 


where e is internal energy, M ini = V a /^yRT 
her, and y is the specific heat ratio. 

The Euler form of the governing equations is easily obtained by 
dropping the second order viscous terms [3] . In this instance, the tan- 
gential and normal velocity components at solid boundaries must be cal- 
culated; this is conveniently handled using simple cosine relations [7] . 


For incompressible flow, the energy equation reduces to the equation 
for the scalar transport of temperature. The pressure is obtained from 
solution of the "discrete" momentum equations and a simple Poisson equa- 
tion based on the SIMPLE algorithm developed by Patankar [8] . A poten- 
tial function, which is solved from the Poisson equation, is used to 
correct the velocity components calculated in the previous time step. 

THE FINITE ELEMENT METHOD 

Bilinear isoparametric quadrilateral elements are used to discretize 
two-dimensional problem domains and trilinear hexahedral elements for 
three-dimensions. The standard weak formulation of the Galerkin weighted 
residual technique is employed to cast Eqs . (1-5) into their integral 

form: 


Conservation of Mass 
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Conservation fif Momentum 
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Conservation of Energy 
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where D denotes the computational domain and W-[_ is the weighting func- 
tion. The p, u, v, and variables are represented by the trial approx- 
imations 
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(16) 


where Nj_ is the linear basis function; in this instance, Wj_ = Nj_. The 
matrix equivalent formulations of Eqs. (11-15) can be expressed as 


A^p + uC 7 p + pC T u = 0 
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(20) 
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M E t + K E t + A (a~) E t = F E 


( 21 ) 


where the * refers to time differentiation, _ denotes column vector, and 
u is the velocity vector. The matrix coefficients are defined (using 
integration by parts for the second order viscous terms) as 
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(26) 
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(29) 

(30) 

(31) 

(32) 


where the i, j, and k subscripts denote summation over the local nodes 
within an element and dS u , dS dS w , dS E represent the boundary segments 
over which u, v, w, and E^ are specified. 

Based on the procedure discussed in Pepper and Humphrey [3], the 
following modifications are employed: 

Mass Lumping 

Equation (22) is modified to the form 
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M = b tJ N t d H (33) 

where 6 LJ is the Kronecker delta. This operation produces a diagonal 
matrix in lieu of the sparse global matrix consisting of n x n nodes, 
thus eliminating the need for implicit matrix solution - hence 

[M]' 1 = 1/Mj. 

Reduced Integration 

The value of an integral over an element is obtained at the element 
centroid using one Gauss point. The shape function becomes 

N i = ^ (2 -0) (34 a) 

A/,-=4 (3-0) (345) 

O 

The advection velocities are evaluated at the centroid of an element 
and multiplied by the average gradient over the element, i.e., 
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The average velocities are obtained as 
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The gradient terms are simple 4x4 matrices in 2-D and 8x8 matrices in 
3-D. For example, c* y is defined for the 2-D quadrilateral as 
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where a = y2 _ Y4 anc * b-y 3 ~y]_. 

The diffusion terms are evaluated in 2-D as 
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and in 3-D 


where 
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The gradient terms have been incorporated to take advantage of the c* Jt 
clj, and cfj formulations. 


Petrov-Galerkin 

The weighting functions associated with the Petrov-Galerkin formula- 
tion are obtained by perturbing the shape functions such that 
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where h e is the element size [3] and a is defined as 


a = coth [3 ( . - 


Pi 


t=l,2 


(43) 


with p j = | u | h e R e /2, which is the cell Reynolds number; for the energy 
equation, p 2 = R r Pi- For the Euler equations, R e nor P r appear in the 
expressions for p. This form of anisotropic balancing diffusion acts in 
the direction of the propagation of the perturbation (velocity field) . 

The precise amount of artificial diffusion (for eliminating the shortest 
waves) and direction in which it must be added for optimizing accuracy 
are calculated for each element. 

Time Integration 

An explicit forward-in-time Euler scheme is used to advance the 
discretized equations in time. For example, u is solved from the rela- 
tion 

u n+1 =u n + &tM~\F u -Ku n - A(u)u n -C x p) (44) 
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Application of this scheme to the continuity, v- and w-momentum and 
energy equations follows similarly. Since the discretized equations are 
solved Explicitly, a stability constraint on the time step must be 
imposed. The Courant and diffusive limits associated with a forward 
Euler scheme are calculated over each element, and the time step 
adjusted to the minimum value within the computational domain. 

EXAMPLES 

A set of six example problems are solved by the FEM to demonstrate 
the versatility of the method. Most of the examples were initially run 
on an IBM PC AT (enhanced with a 33MHz accelerator board [9]) and an 
Everex 25MHz 386 machine; the first four problems required nodal meshes 
using less than 2000 nodes. Although problems requiring upwards of 5000 
nodes have been successfully run on PCs [2], such jobs may require many 
hours of calculational effort. 

The example problems were also run on an Alliant VFX/40 mini- 
supercomputer. A Silicon Graphics Personal IRIS workstation served as a 
graphical front end to the Alliant. This combination of hardware allowed 
near real time graphical displays of the transient solutions as they 
converged to steady state. Convergence was assumed to be reached when 
differences between old and new values (all variables) equalled lO - ^. 

Example 1 : 

The first problem deals with fluid recirculation within a square 
cavity (0 < x < 1 ; 0 < y < 1 ), as shown in Fig. l(a,b). The top surface repre- 
sents a fluid moving horizontally left to right at u=l, v=0 with R e = 
5000; the remaining three walls are no-slip surfaces. 

A non-uniform rectangular mesh of 51x51 nodes is used to model the 
fluid motion. Ghia, et al [10] and Gresho, et al [4] also analyzed flow- 
in a cavity for a large range of Reynolds numbers. Figure 1(b) shows the 
streamline distributions, and agrees with their published results. A 
primary vortex develops within the center region of the cavity; second- 
ary vorticies occur at the two bottom corners and near the upper left 
corner . 

Example 2 : 

In the second problem, natural convection within a rectangular 
enclosure is simulated for R a =10^ and P r =l . A 20x20 nodal mesh is used 
with zero initial conditions for u, v, and T (Fig. 2(a)). 

At R a =105, two distinct cells develop from one large single cell 
shortly after the calculation is begun. This recirculation pattern is 
consistent with multicellular development as the Rayleigh number 
increases. Fig. 2(b,c) shows the streamlines and isotherms for steady 
state conditions, and agree with results in the literature [11], At 
higher Rayleigh numbers, a finer grid is required near the walls to 
properly account for the boundary layer growth. 

Example 3 : 

The ability of the FEM to model convective cooling of a heated block 
is demonstrated in the third example problem. Prediction of the recircu- 
lation zone behind the block and the thermal plume eminating from the 
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block are examined for R e =P e =2000. The problem geometry and mesh are 
shown in Fig. 3(a,b); the mesh consists of 732 nodes (665 2-D elements). 

Since the mesh is coarse, spurious results are expected for 
R e -P e =2000; at lower values, the mesh is adequate. In this case, the 
ability of the Petrov-Galerkin scheme to control dispersion errors (os- 
cillations in the solution) is assessed. Figure 4(a,b,c,d) shows veloc- 
ity vectors and isotherms (with and without Petrov-Galerkin weighting) . 
At high Reynolds numbers, the FEM yields poor velocity distributions 
without Petrov-Galerkin weighting; the temperature distribution shows 
appreciable dispersion - such solutions generally indicate inadequate 
mesh refinement. Use of Petrov-Galerkin reduces the oscillations and 
provides a good solution with a less than optimal mesh. 

Example 4 : 

In the fourth problem, a planar 15° ramp sits on the lower wall of a 
two-dimensional duct. The inlet Mach number is 2.28. Based on 1-D ideal 
analysis (weak solution for attached shock waves) , an oblique shock 
develops at an angle of 40° relative to the freestream. Figure 5(a,b,c) 
shows the flow schematic and problem geometry. The shock reflects off 
the top wall at an angle of 40.5°. In region 2, M 2 = 1.69 and P 2 /P 1 
(static) = 2.34. In region 3, M 3 = 1.15 and P 3 /P 1 = 4.91. 

Figure 6 (a,b) shows density contours and velocity vectors for the 
viscous solution. A time step of A(=10 ' 4 is initially required. The 
boundary layer development is apparent; at the upper wall, the reflected 
shock interacts with the boundary layer and creates a small recircula- 
tion zone. For this viscous solution, R e = 10^; turbulence has not been 
included in the model as yet. Although the shock is smeared, the shock 
angle and downstream Mach numbers agree with the 1-D analysis. A more 
optimized grid refinement (more uniformity and nodes) is required to 
reduce spreading of the shock. Using 2x2 Gauss points versus one Gauss 
point increases accuracy, but also increases the compute time. 

Example 5 : 

The fifth problem is concerned with supersonic flow over a projec- 
tile within a two-dimensional chamber. This problem is similar to the 
Scramaccelerator work conducted by Pratt, et al [12] and Humphrey [13]; 
the projectile is designed so that the reflected shock off the tube wall 
initiates combustion (oblique detonation wave) , thereby producing an 
increase in pressure behind the shocks and a net thrust forward. In this 
example, the upper wall moves at u=l (non-dimensional) . Approach flows 
for Mach 2.68 (nose angle 20°) and Mach 5 (nose angle 25°) are simulated 
with the initial flow field assumed to be uniform (i.e., equal to the 
inlet conditions) throughout the domain. The problem geometry and mesh 
are shown in Fig. 7. 

Figure 8 (a,b,c,d) shows pressure and Mach contours for Mach 2.68 and 
Mach 5 inlet conditions, respectively. An oblique bow shock forms over 
the nose and reflects off the top wall. High pressure gradients develop 
in the region between the wall and the projectile. Figure 9(a,b) shows 
an expanded view of the velocity vectors in front of the projectile. In 
Fig. 9(a), boundary layer separation occurs on the upper part of the 
nose as a result of the interaction with the reflected shock; in this 
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case, the inlet flow speed is not high enough and the shock inhibits 
forward motion. At Mach 5 and a 25° leading edge (Fig. 9b), the oblique 
shock reflects off the upper wall and strikes the projectile just behind 
the shoulder (near optimal performance) . A weak shock train is created 
in the thin region between the upper wall and projectile body; boundary 
layer separation occurs along the top of the projectile body. The pre- 
dicted oblique shock angles and downstream flow conditions near the 
front shoulder of the projectile agree quantitatively with 2-D results 
obtained by Pratt, et al [12] using a modified SALE algorithm [14] . 

Example 6 : 

For the last problem, 3-D flow is calculated around a set of heated 
obstacles which are convectively cooled by cold air. The physical domain 
and mesh are shown in Fig. 10 (a,b) . The mesh consists of 2868 hexahedral 
elements (which easily permits execution on a Personal IRIS worksta- 
tion) ; the Reynolds number is R e =lC)3 and P r =l . 0 . This type of problem is 
commonly encountered in the computer industry where cooling of computer 
chips is critical. 

Figure ll(a,b) gives normal and perspective views of the 3-D veloc- 
ity vectors within the channel. Recirculation of the flow occurs behind 
the blocks, and small secondary cells develop in the corners. The 
pressure and isotherm distributions are shown in Fig. 12 (a, b) for the 
x-y plane near the top of the channel. Distinct thermal plumes eminate 
from the heated blocks; plume impingement from the left forward block 
occurs on the small mid-stream block. It is well known that when flow 
separates at the corners of blocks, horseshoe-like vortices are gener- 
ated [15] . Further work is underway to examine this phenomena using a 
28,000 hexahedral element mesh on the Alliant. 

CONCLUSIONS 

A finite element method using simple modifications is used to solve 
2- and 3-D problems for compressible and incompressible flows with heat 
transfer. The modified FEM employs equal order basis functions for all 
unknown variables, mass lumping, reduced quadrature, and Petrov-Galerkin 
weighting; transient solutions are solved using an explicit Euler 
scheme. Formulation of the local (and global) matrices are simple, and 
solution speeds are quick. The technique yields accurate results for a 
wide variety of flows providing the elements are not too distorted. The 
method appears attractive for workstation class machines. 

ACKNOWLEDGEMENTS 

We wish to thank Mr. Andy Singer, Advanced Projects Research, Inc., 
for his help in generating many of the graphical displays, and Dr. Juan 
Heinrich' and Mr. Frank Brueckner, University of Arizona, Tucson, Ari- 
zona, for their assistance in generating the Scramaccelerator results. 


321 



REFERENCES 


1. Thompson, J. F., "Numerical Solution of Flow Problems using Body 

Fitted Coordinate Systems," Computational Fluid Dynamics , 1, (W. Koll- 

man, ed.) , Hemisphere Pub. Co., Washington, D.C. (1980). 

2. Pepper, D. W. and Singer, A. P., "Calculation of Convective Flow on 
the Personal Computer using a Modified Finite Element Method, " to appear 
in Num. Heat Transfer (1990) . 

3. Pepper, D. W. and Humphrey, J. W., "A Hybrid Finite Element Method 
for Compressible Flow," AIAA Paper 90-0399 (1990). 

4. Gresho, P. M., Chan, S. T., Lee, R. L., and Upson, C. D., "Modified 
Finite Element Method for Solving the Time-Dependent Incompressible 
Navier-Stokes Equations, Part 1: Theory," Int . J. Num. Methods in Flu- 
ids, 4, 557-598 (1984) . 

5. Yu, C. C. and Heinrich, J. C., "Petrov-Galerkin Methods for the Time- 
Dependent Convective Transport Equation," Int. J. Num. Methods in Eng., 
23, 883-901 (1986) . 

6. Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., Computational 
Fluid Mechanics and Heat Transfer, Hemisphere Pub. Co., Washington, D. 

C. (1984) . 

7. Shapiro, R. A., "An Adaptive Finite Element Solution Algorithm for 
the Euler Equations, " MIT Report CFDL-TR-88-7 (1988) . 

8. Patankar, S. V. , Numerical Heat Transfer and Fluid Flow, Hemisphere 
Pub. Co., Washington, D.C. (1980).- 

9. Definicon Systems, Inc., "PM Series Product Manual," Definicon Sys- 
tems, Inc., 1100 Business Center Circle #5, Newbury Park, CA 91320 
(1989) . 

10. Ghia, U., Ghia, K., and Shin, C., "High-Re Solutions for Incompress- 
ible Flow Using the Navier-Stokes Equations and a Multi-grid Method, " J. 
Comp. Phys., 48, 387-401 (1982). 

11. Jones, I. P. and Thompson, C. P. (ed.), "Numerical Solutions for a 
Comparison Problem on Natural Convection in an Enclosed Cavity, " AERE-R 
9955, Harwell, England (1980) . 

12. Pratt, D. T., Humphrey, J. W., and Glenn, D. E., "Morphology of a 
Standing Oblique Detonation Wave," AIAA Paper 87-1785 (1987). 

13. Humphrey, J. W., "Study of an Oblique Detonation Wave Ramaccelerator 
Driven Hypersonic Test Facility," APRI Final Report, NASA Langley Con- 
tract NAS 1-1 8 8 02 (1989) . 

14. Amsden, A. A., Ruppel, H. M., and Hirt, C. W., "SALE: A Simplified 
Arbitrary Lagrangian-Eulerian Computer Program for Fluid Flow at All 
Speeds," Los Alamos Nat. Lab. Report LA-8095, UC-32 (1980). 

15. Hunt, J. C., Abell, C. J., Peterka, J. A., and Woo, H., "Kinematical 
Studies of the Flows around Free or Surface Mounted Obstacles; Applying 
Topology to Flow Visualization," J. Fluid Mech., 86, 1, 179-200 (1978). 


322 



u = 1 , \/=0 



u = v =0 


a.) Mesh 



b. ) Streamlines 

Figure 1= Flow in a Cavity - Re=50G0 
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b. ) Mesh 


Figure 3 0 Forced Convective Cooling - Re=Pe=200Q 
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a. ) w/o Petrov-Galerkin Weighting 



b. ) Temperature Distribution 
(w/o Petrov-Galerkin Weighting) 


Figure 4. Forced Convective Cooli 
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c. ) Velocity Vectors 
(w/ Petrov-Galerkin Weighting) 



d. ) Temperature Distribution 
(w/ Petrov-Galerkin Weighting 


- Velocity Vector and Isotherms 




u ■= v = <3T / <3/7 «* 0 


p-1 outflow — 



c. ) Mesh for viscous solution 

Figure 5. Supersonic flow over a 15° Ramp 
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Moo = 2.68 



Moo = 5.00 



a. ) Pressure Contours 


Moo = 2.68 




b. ) Mach Contours 


Figure 8. Scramacceierator 
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Moo = 5.00 


Figure 9. Inlet Velocity Distribution - Scramaccelerator 
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mi min 


u = 1 


v = 0 


w = 0 


T = 0 




b. ) Mesh 

Figure 10. Channel with Obstacles - Forced convective Flow 
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a. ) View from Above 



b. ) Perspective View 

Figure 1 1 . Velocity Distribution 
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a. ) Pressure 



b. ) Isotherms 

Figure 12. distribution in X-Y Plane 
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